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Abstract 

Stochastic solutions are obtained for the Maxwell- Vlasov equation 
in the approximation where magnetic field fluctuations are neglected 
and the electrostatic potential is used to compute the electric field. 
This is a reasonable approximation for plasmas in a strong external 
magnetic field. Both Fourier and configuration space solutions are 
constructed. 

PACS: 52.25Dg, 02.50Ey 
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1 Introduction. The notion of stochastic so- 
lution 

The solutions of linear elliptic and parabolic equations, with Cauchy or 
Dirichlet boundary conditions, have a probabilistic interpretation. These are 
classical results which may be traced back to the work of Courant, Friedrichs 
and Lewy [1] in the 1920's and became a standard tool in potential theory[2] 
[3] A simple example is provided by the heat equation 

1 ^2 

5iM(t,x) = -— With u{0,x) = f{x) (1) 
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with solution written either as 

« (t, x) = -i= I i= exp (^-^^^ j / (y) dy (2) 

or as 

uit,x)=EJiXt) (3) 
denoting the expectation value, starting from x, of the process 

dXt = dWt 

Wt being the Wiener process. 

Eq. ([1]) is a specification of a problem whereas ([2]) and ([3]) are solutions in 
the sense that they both provide algorithmic means to construct a function 
satisfying the specification. An important condition for ([2]) and ([3]) to be 
considered as solutions is the fact that the algorithms are independent of 
the particular solution, in the first case an integration procedure and in the 
second the simulation of a solution-independent process. In both cases the 
algorithm is the same for all initial conditions. This should be contrasted 
with stochastic processes constructed from particular solutions, as has been 
done for example for the Boltzman equation 

In contrast with the linear problems, explicit solutions for nonlinear par- 
tial differential equations, in terms of elementary functions or integrals, are 
only known in very particular cases. Hence it is in the field of nonlinear 
partial differential equations that the stochastic method might be most use- 
ful. Whenever a solution- independent stochastic process is found that, for 
arbitrary initial conditions, generates the solution in the sense of Eq. ^ , an 
exact stochastic solution is obtained. In this way the set of equations for 
which exact solutions are known might be considerably extended. 

The exit measures provided by diffusion plus branching processes [6] [7] 
[8] [9] [To] [11] as well as the stochastic representations recently constructed 
for the Navier-Stokes[I2j ^ISj ^ [15] [IB] [I?], the Poisson-VlasovfTg] [12], 
the Euler[2Dj and a nonlinear fractional differential equation[21J define initial 
condition-independent processes for which the mean values of some function- 
als are solutions to these equations. Therefore, they are exact stochastic 
solutions. 

Typically, in the stochastic solutions, one deals with a process that starts 
from the point where the solution is to be found, the solution being then 
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obtained from a functional computed along the whole sample path or until 
the process hits the boundary. In addition to providing new exact results, 
the stochastic solutions are also a promising tool for numerical implementa- 
tion. Here the relevant question is to know when a stochastic algorithm is 
competitive with the existing deterministic algorithms. Although there is no 
general answer to this question, there are a few considerations that suggest 
where and when stochastic algorithms might be useful, namely: 

(i) Deterministic algorithms grow exponentially with the dimension d of 
the space, roughly N'^ (-^ being the linear size of the grid). This implies 
that to have reasonable computing times, the number of grid points may not 
be sufficient to obtain a good local resolution for the solution. In contrast a 
stochastic simulation only grows with the dimension of the process, typically 
of order d. 

(ii) In general, deterministic algorithms aim at obtaining the behavior of 
the solution in the whole domain. That means that, even if an efficient deter- 
ministic algorithm exists for the problem, a stochastic algorithm might still 
be competitive if only localized values of the solution are desired. This comes 
from the very nature of the stochastic representation processes that always 
starts from a definite point in the domain. According to what is desired, 
configuration or Fourier space representations should be used. For example 
by studying only a few high Fourier modes one may obtain information on 
the small scale fluctuations that only a very fine grid would provide in a 
deterministic algorithm. 

(iii) Each time a sample path of the process is implemented, it is indepen- 
dent from any other sample paths that are used to obtain the expectation 
value. Likewise, paths starting from different points are independent from 
each other. Therefore the stochastic algorithms are a natural choice for par- 
allel and distributed implementation. Provided some differentiability condi- 
tions are satisfied, the process also handles equally well simple or complex 
boundary conditions. 

(iv) Stochastic algorithms may also be used for domain decomposition 
purposes [22] [23] [21]. One may, for example, decompose the space in sub- 
domains and then use in each one a deterministic algorithm with Dirichlet 
boundary conditions, the values on the boundaries being determined by a 
stochastic algorithm. 

Stochastic solutions also provide an intuitive characterization of the phys- 
ical phenomena, relating nonlinear interactions to cascading processes. By 
the study of exit times from a domain they also provide access to quantities 
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that cannot be obtained by perturbative methods [25] [SB] . 

One way to construct stochastic solutions is based on a probabihstic in- 
terpretation of the Picard series. First, the differential equation is written as 
an integral equation. Then there are two possibilities. 

In the first case the series is rearranged in a such a way that the coef- 
ficients of the successive terms in the Picard iteration, including the initial 
condition term, obey a normalization condition. The stochastic solution is 
then equivalent to importance sampling of the normalized Picard series. In 
this case the stochastic process, that constructs the solution, is a branching 
process, the branching being controlled by the nonlinear part of the equation. 

The second possibility occurs when the initial condition term is not mul- 
tiplied by a probability factor. Then, the integral of the integral equation 
may still be given a probabilistic interpretation but the process that is used 
for the construction of the solution is more general tree-indexed stochastic 
process. 

In this paper, pursuing the work on kinetic equations initiated in [18] and 
[19j . solutions are obtained for the Maxwell- Vlasov equation in the approxi- 
mation where magnetic field fluctuations are neglected and the electrostatic 
potential is used to compute the electric field. This is a reasonable approx- 
imation for plasmas in a strong external magnetic field. Both Fourier and 
configuration space solutions are constructed. 

In Sect. 2.1 one discusses the formulation of the full Maxwell- Vlasov sys- 
tem as an integral equation for the charge densities. In Sects. 2.1 to 2.4 the 
solutions to the Fourier-transformed equation are obtained both for a static 
uniform and a slowly varying magnetic field. The stochastic processes asso- 
ciated to the construction of these solutions are branching processes with the 
waiting time controlled by the velocity Fourier component and the densities 
(anti-) evolved by the linear part of the equation. 

In Sect. 3 one deals with the configuration space equation and in this 
case the most natural stochastic formulation involves a general tree-indexed 
stochastic process. 
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2 The Poisson-Vlasov equation in a magnetic 
field 



2.1 The Maxwell- Vlasov system 

Consider a two-species Maxwell- Vlasov system in 3-1-1 space-time dimensions 

^+ V -VJ, + ^(e+-xb]- VJ, = (4) 
at rrii \ c I 

{i = 1,2), with 

E-cV^xB = -in^ei / v fidh 



dt 

d_ 
dt 



— B+cV^xE = 



V,-S = (5) 

To study the Vlasov- Maxwell system as a nonlinear equation for the /j {t, x, v) 
densities one has to obtain explicit expressions for the electromagnetic fields 
in terms of the charge densities. Define scalar and vector potentials 

c at 

B ^ VxA (6) 
which, in the Lorentz gauge (V- A = 0), obey the equations 

i 

^A--,^ = vM'v (7) 

i 

Using the (retarded) Green's function for the wave equation. 



one obtains 



A {x,t 



3^/ 



X — X 




dX9) 



In writing ([9]) as a source term solution of ([7]), one has assumed that the initial 
conditions for the equations ([7j) vanish together with their time derivatives 
or, alternatively, that the initial time is in the remote past so that there 
are no more contributions from the initial conditions. For a more general 
solution which should be used for transitory phenomena or plasma probing 
by short localized pulses see [27] . 
Use of ([6]) yields 



E 



+ 



c| X — X I 



VV X 



+ 



dtj f \^x ,v,t 
TTjdt \ f (x',v,t- L- 



(10) 

Replacing now Eq. (|TOl) in (jlj) one sees that, as a function of the densities 
f{x,v,t), Maxwell- Vlasov is a nonlinear and nonlocal (in time) differential 
equation. Its full stochastic solution treatment will be dealt with elsewhere. 
Here one deals with approximations of practical interest for fusion plasmas. 
First notice that for non-relativistic plasmas the last terms in ( fTOl) are small 
and in the quasi-static approximation one has 



E = J d^x'^^Ci ^ 3 / d^vf{x',v,t) 



B = I d'x' Er-r^ I 



V \ V X X 



a;' I I / {x',v,t) 



Furthermore, for microturbulence studies in fusion plasmas in strong mag- 
netic fields, a reasonable approximation neglects magnetic field fluctuations 
and uses the electrostatic potential of the charges to compute the electric 
field. This is what will be called Poisson- Vlasov in a static (external) mag- 
netic field. 
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2.2 Poisson-Vlasov in a static magnetic field 



In this approximation the equation is 

= ^+f^-V. + — -xi(x)-V.l/. 
at \ rrii c 

+ — /rfVVej d\fj (x\u,t) — ^ - VJ^ix^v^t) (12) 
"^i J J ^ ^ \x — X \ 

or, in the Fourier transformed version 

^(e,t) = ^/ A/(^,t)e^^"' (13) 

with 7/ = {x, and ^ = (^2, = (^i, ^2), 



dt 



(14) 



The aim is to obtain stochastic solutions both for the equation in configu- 
ration space and for its Fourier-transformed version. Because of the locahzed 
nature of the stochastic solutions, as discussed in the introduction, both types 
of solutions are useful for the applications. If, in a plasma confinement de- 
vice, one is interested in the behavior of the solution at a particular point 
(for example at a point either in the core or in the scrape-off layer) then it is 
the solution in configuration space that is useful. If however one is interested 
in the nature of the turbulent fluctuations it is probably the study of high 
Fourier modes in the Fourier-transformed equation that will be most useful. 

Eqs.( ll2p and fll4p have linear and nonlinear parts. The linear evolutions 
are, respectively 

(x,t;,t) = e-^Q-'/f (x,t;,0) 

the operators Qr, and being 

Qn = V ■\/^ + X B (x) -V^ 

Q( = -ri-V6-^V5,x5(-zVgj-f2 
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(16) 



From fllSI) and 0161) it follows that the linear evolution of the function argu- 
ments X, V, ^1 and ^2 is ruled by the following equations 



(17) 




^1 (t) + ^i?HV^, W)x6W 



(18) 



and 



iVa (t) 
|V,, (t) 



(t) 

^(Ve, (t)xi?(-zVa (t))) 



(19) 



One sees from flTTl) that the linear evolution of the densities f {x,v,t) 
in configuration space acts only on the arguments of the function. How- 
ever, from (fT8|) and (fT9|) one also sees that, if the magnetic field B is not 
constant in space, the linear evolution of the Fourier transformed densities 
^i^^ (^1,^2,^) is more complex, involving derivatives of the Fourier density. 
For the stochastic solutions, one associates a process to each function, there- 
fore it is not convenient to use the full linear part of the evolution operator in 
the Fourier-transformed equation. Such problem does not exist if the static 
magnetic field is also uniform in space. Here one starts by studying this case, 
which is then extended to the case of a slowly varying magnetic field. 

2.3 Fourier-transformed Poisson-Vlasov in a static uni- 
form magnetic field 

Consider a uniform magnetic field B=Bo= Bq el. In this case it is possible 
to obtain an explicit form for the evolution of the linear part. 




) 



(20) 



that is 




^1 (^) + ^5ox6 (t) 







(21) 
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with solution {t) =^i and 
(6 W)^ = - {^2) y sin Uit + {^2)^ cos Uit + |(^i)^sinwit + (^i)^(cosw,t- 1)| 

(22) 

with cij,- = ii^Q inverse relation beine 

' crrii ' 

fe), \ _ ( coscj.t sincuit \ / (6 (^))z - i ((6).sin^,t + (^1)^ (cos^^t - 1)^ 
" ^-sin^.t cosa;.t J (^^(i))^_i.|(^^)jl_eosa;,t) + (ei),sin^/ 

(23) 

In integral form Eg. (1141) becomes 



(6,6, t) = F, (6,6 (t) , 0) - ^ / rfs / (6 - el, 6 (s) , t - s 

mi Jo J \ 

-^§^j:hF,{('..0,t-s) (24) 



A stochastic solution is going to be written for the following function 

X.(ei,6,t) = e-'-(l«^l)^^i|liM (25) 

where 7(|6l) = 1 if |6l ^ 1 and 7(|6l) = 1^21 otherwise. /i(6) a positive 
function to be specified later on. The integral equation for Xi {^i,^2,t) is 

rui h (6 ) 

* 7 (16(g) I) Jt-shi\i2(s)\)-tM\i2\) f J3 



X / ds—— — - — -e 

N{^i,^2,t) 



A 
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with \Q , 



and 





-1 


h (6 - Q h (el) 


( 




^h*h 





(27) 



(28) 



Eq. (|26|) has a stochastic interpretation as an exponential plus a branching 
process. The survival probability up to time t of the exponential process is 



g-<7(l6l) 

and dsU {^i,C,2,s) is the decay probability in time ds, with 

^_ 7(|6(^)|)e(*-)^(fe(-)l)-*^(fel) 

(^1, ^2; t) being a normalizing function 

1 



(29) 



(30) 



1 



-47(161) 



ds-f (IG (s)|) e(*-^)^(l«2(^)l)-*^(l«2l) (31) 



In the branching process, p li^^i is the probability that, from a 

mode, one obtains a (^i — ^i, ^i) branching with ^[ in the volume [^[^ ^'i + d^^'i) . 

The stochastic interpretation of Eq. (l26|) provides a way to compute the 
solution. Xi (^1; ^2, t) is computed from the expectation value of a multiplica- 
tive functional associated to the process. Convergence of the multiplicative 
functional hinges on the fulfilling of the following conditions : 



(A) ™ < 1 

(B) (|e;r'A*A)({i) </!({,) 

Condition (B) is satisfied, for example, for 



(1 + 161^)' 



and 



c < 



3n 



(32) 
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Indeed computing ^ h* h one obtains 



2|?l|'(l6l'+4 



(33) 
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Then ^|^^| ^ /i * /ij (^i) is bounded by a constant for all |^i|, and choos- 
ing c sufficiently small, condition (B) is satisfied. 

Once h{^i) consistent with (B) is found, condition (A) only puts restric- 
tions on the initial conditions. Now one constructs the following backwards- 
in-time process, denoted X {^i,C,2,t): 

Starting at (,^1, t), a particle of species i lives for a 11 (^1, ,^2, s) —distributed 
time s, up to time t — s, with survival and decay probabilities given by fl29l) 
and (150]) . At its death a coin 1^ (probabilities |, is tossed. If = two new 
particles of the same species as the original one are born at time t ~ s with 
Fourier modes (^1 — ^2 (s)) and [^i, O) with probability density p (^1, ^i) . 
If Is = 1 the two new particles are of different species. Each one of the new- 
born particles continues its backward-in-time evolution, following the same 
decay and branching laws. When one of the particles of this tree reaches 
time zero it samples the initial condition. The multiplicative functional of 
the process is the product of the following contributions: 

- At each branching point where two particles are born, the coupling 
constant is 

^ y \ _ 8vre.e,Ar(ei,6,t) {l^'i^ ^ * h) (^i) g (,) . 

- When one particle reaches time zero and samples the initial condition 
the coupling is 

The multiplicative functional is the product of all these couplings for 
each realization of the process X(^i,^2,^)- The solution Xi(^i5^2,^) is the 
expectation value of the multiplicative functional. 

Xi (6, 6, t) = E |n (^gogQ ■ ■ ■ j (gug'.i ■ ■ ■ ) (^gtjg'ij ■ ■ ■ ) } (36) 
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Fig.l illustrates a realization of the process. Notice that the label (53 — Si) 
denotes the mode ^2 = ( ant i-) evolved during the time S3 — s-[ according to 
Eq.dl] 



time 



1 1 3 r' 



■1 ■■i'--2' 1 



o 



6 6 



o Fourier 
mode 



Figure 1: A sample path of the stochastic process X (^1,^2,^) 
The process itself is the limit of the following iterative process, 



= X^{^l,^2{t),0)^>r]+9n[^l,^'l,s'^ (37) 

xXf ) (6 - C 6 {s),t- s) Xf ) (e;, 0, t - s) l[,.=o] 

+^7^, (^1, el, ^.^'^ (6 - ^'1, 6 (s) , t - s) x]'^) (e;, 0, t - s) i[.<.]i[^,=i] 

with the solution being 

X. (6,6,t)=E{iimxf 
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With the conditions (A) and (B) and choosing the c constant in h (^i) 
such that 



mini {mi} h (^i 



< 1 (38) 



the absolute value of all coupling constants is bounded by one. The branch- 
ing process, being identical to a Galton- Watson process, terminates with 
probability one and the number of inputs to the functional is finite (with 
probability one). With the bounds on the coupling constants, the multi- 
plicative functional is bounded by one in absolute value almost surely. 

Once a stochastic solution is obtained for Xi {^i,^2,t), one also has, by 
fl24|) . a stochastic solution for Fj (,^1,^2,^)- Summarizing: 

Theorem 1 The stochastic process X {C,i,^2,t), above described, provides 
through the multiplicative functional a stochastic solution of the Fourier- 
transformed Poisson- Vlasov equation in a uniform magnetic field for arbi- 
trary finite values of the arguments, provided the initial conditions at time 
zero satisfy the boundedness conditions (A). 

2.4 Fourier-transformed Poisson- Vlasov in a static non- 
uniform magnetic field 

The result is now generalized to the case of a static non-uniform magnetic 
field. Decompose the Fourier transform of the magnetic field into 

B{ii) = {2'nf'Bo5\ii)+l>{ii) (39) 

where Bo might be the average of the field in a region of interest and the 
non- uniform part h (G) is assumed to be small, in a sense to be specified 
later. Then the integral equation becomes 

F, (^1, 6, t) = F, (ei, 6 it) ,0)-^ fds [ d%F, - ^2{s),t-s 

mi Jo J V 

X 6 (s) ■ p (Q X V6{.)) (6 - is),t- s) (40) 
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where, as before, the dynamics of the arguments ^2 (t) in Fj is controlled by 

the constant Bo (Eqs. lpT]) and (!22|) ). A stochastic solution will be obtained 
for the function 

X.(ei,6,t) = e-*-(l«^l)^^%%^ (41) 
with integral equation 



xe(*-'*)'>(l€2WI)-*7{16l) 



j|x.(ei-eU2(s),t-s) (42) 



As before, a backwards- in-time process, rooted at {^1,^2, t), is considered. 
The survival and branching probabilities are also ruled by (!29|l and (130|) . 
However now, whenever the propagating particle dies, there are three distinct 
possibilities. Either two new particles of the same species (or of opposite 
species) are born at time t — s with Fourier modes (^1 — ^2 (s)) and (^;^, O) 
with probability density p (^1,^1) given by (l28!l . or it is just one particle with 

mode (^1 — ^2 (s)) that is born and the process samples the field b {^i}- 
That is, the particle samples the non-uniform field and is scattered by it. 
This particle also receives an operator label 

The operator labels are subsequently inherited by the offspring of this parti- 
cle and accumulate until they are finally applied to those offspring particles 
that reach time zero. There is no ambiguity in the application of the op- 
erators to the final particles because both the ^2 [t) argument of the final 
particles and the derivatives V6(s) should be expressed in terms of the ini- 
tial (^1, ^2, VCi) V6 using the solutions of Eqs. fllSHlQl) . 

Denote by Y {C,i,C,2,t) the process obtained as the iterative limit of the 
construction described above. A realization of the process is illustrated in 
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Fig. 2. The boxed K ('£",£5 fs^)) denotes the operator label that is attached 
to thi 



time 



",0 ) 



(t-^"'0(s -s )) 
1 1 3 



'1 ^i'^2^ 1 



O 



O 



o 



Fourier 

mode 



Figure 2: A sample path of the stochastic process Y (^1, t) 

The solution Xi {^i^^2,t) of the equation is then obtained from the aver- 
age value of a multiplicative functional associated to the process. For each 
realization of the process, the functional is the product of the following fac- 
tors: 

- At each branching point where two particles are born , the coupling 
constant is 

^ , _ 167re.e,iV (6,6, t) (kH'' ^ * ^) fe) g (.) ■ ^[ 

m M^i) 7(l6(.)|) ^ ^ 

When one particle reaches time zero and samples the initial condition 
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the coupling is 

90^ (^1, 6) = K 6 K [Ci2 {S.)) ■■■K (er, 6 {Sn)) ^^^f^ 

' (45) 

In addition to condition (B) of the previous subsection, a sufficient con- 
dition for the convergence of the functional is 



mini {m} 



< 1 



and 



(46) 



< 1 (47) 



for arbitrary n and arbitrary values of the arguments ^i,^2{s)- This last 
condition requires boundedness and smoothness of the initial condition as 

well as a sufficiently small (as compared to h{C,i)) non- uniformity field b (^i)- 
Summarizing: 

Theorem 2 The stochastic process Y {^i,^2,t), above described, provides a 
stochastic solution to the Fourier-transformed Poisson-Vlasov equation in a 
static non-uniform magnetic field, provided the initial conditions at time zero 
and the non-uniform part of the field satisfy the condition ( [77| j. 



2.5 The configuration space equation 

Eq.f ll2l) in integral form is 



ft {x,v,t) 



e-'^''f^{x,v,0)- ^ dse-'^"^ j d^x^Cj j d\fj(^x',u,t 



or 



\x — X 



fi{x,v,t) = f,{x{t),v{t),0)- 



nii 



[ ds J d^x Cj J d^ufj (^x ,u,t 



— s 



^ rrr — ^ ■ ^Hs)fi {x (s) ,v{s),t-s) 

\x (S) — X \ 



(49) 
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x{t), V (t) being the solutions of (fTTI) with initial conditions. 

In the Fourier-transformed equation, division by 7(1^2!) not only regu- 
larizes the velocity gradient as it also, through multiplication of Fi (^1,^2,^) 
by e~*'''*^'^^'\ introduces a natural time scale for the exponential process that 
controls the branching. Here, because division by Vt, does not make sense, 
there is no natural exponential time scale. One could nevertheless multiply 
fi (x, V, t) by e"^*, with A a constant, as in Ref. [2U] for the equation without 
magnetic field. However, because of the nonlinear nature of the second term 
in poj) . this introduces strong limitations on the range of t for which the 
solution may be constructed. Here a different procedure will be followed. 
The price to pay is that, instead of a simple branching process, one needs a 
more complex tree-indexed stochastic process. 
Let Gi ( x,v,t] be the function 

fi (x, V, t) 

Gdx,v,t) = ^^ (50) 



(Piix {t),v (t) 



the ip (^x, f j's being functions to be specified later and (^x (t) , v (t)j the 
function arguments (anti-)evolved by f[T7|l . One obtains the following integral 
equation for Gj ^x, f , 



Gi I X , V 



^ j "'O 

X / d?x'd?up"^^y^^ ( x', M, s j ( x', u, t — s j ^x (s) — x'j 



(pi ( X {t),v (t) 



■^v{s)^i {x (t) , V (t) j Gi (x {s) ,v{s),t 



with V = ri and 

^ \y\ 



(j) 

Px,v,t I X , U, S 



(fj I x' {t — s) , u (t — s) 



A 



x,v,t 



X (s) 



X 



(52) 
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a probability in the space [0,t] x x M^, A^^^^t being the normahzation 
constant 



ft r r [x' (t-s), u {t - s) 

Ai'lt = ds d?x'd'u—^ — ^ (53) 

■^0 x{s)-x' 

One of the simplest choices for the functions ^x, f ^ would be to make it 
proportional to the initial condition 

{x, = kfi V, j (54) 

Then, the probabilistic interpretation would require finiteness of 







t . kfj {x'{t-s),u{t-s),0 

Axl,t,s = I ds I d^x'd^u ■ — ^ (55) 

X (s) — x' 



a quantity that has the nature of a retarded field intensity generated by 
the initial condition. However, the general result will be stated without 

committing to a particular choice oi {^x ^ . 

From Eq.( |5T|) one sees that because the term Gi {^x (t) , v (t) , is not 
multiplied by a probability factor one cannot simply interpret the construc- 
tion of (^x,v, as importance sampling of the Picard series. Nevertheless, 
a probabilistic interpretation may be given through the following tree- indexed 
stochastic process Z (^x, v,t^: 

Rooted at (x,u,t], a, particle of species i propagates backwards- in-time 



until a time t — s when, controlled by the probability p^J\ ^ [ x',u,s ] it gives 



birth to two new particles. One of them is of the same species i and the other 
of the same or the opposite species with probability |. The first particle has 



coordinates (^x (s) , v (s)j and the other coordinates yx', uj determined by 
the probability p^lt { x',u,s]. The first particle also receives an operator 
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label 



K (s) = (x {s) - x') • ^ ^ ^V,(,)<^, (x (t) , V (t)) (56) 

to be subsequently applied to all of its offspring. The original particle, the 
one that gave birth to the two new ones, does not die and proceeds its free 
propagation until time zero. Then, each one of the newly created particles 
has an evolution analogous to the progenitor and during the its evolution 
the operator labels that they inherit at the birth of each new pair are ac- 
cumulated, until they are finally applied to the initial condition when each 
one of the particles reaches time zero. A realization of the process is illus- 
trated in Fig.3. The flags, denoted si, S2, ■ ■ ■ , stand for the operator labels 




The main differences from the Fourier-transformed case are: 
- The progenitor particles never die. 
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- The solution of the equation is obtained from the average over reahza- 
tions of the following quantities 



Gi [x,v,t 



(57) 



2 6 -J 



mi 



A%,tGj [x\u,t-s\K{s) G^ [x{s),v {s),t 



with Gj \^x', u,t — s j and Gi (s) ,v (s) ,t — sj computed in the same 

way until t — s reaches time zero. For each realization the process runs 

from time t to zero. However, the calculation of the quantities Gi for each 
realization runs the opposite way, from time zero to time t. Qualitatively, 
what the process does is to replace the calculation of the integrals in fl5T|) by 
the generation of a family of probability measures and each value of fl57|) is 
a sampling of the corresponding Picard iteration. 

Assume that, with probability one, the iteration ( l57Il converges for all 
realizations of the process. Then the solution of (1511) is obtained from 

, v,t^=E jCi (x, V, } (58) 

Hence, existence of the stochastic solution depends on the boundedness 
and convergence of the iteration in (157|) . Let 



Gi [x,v,0 



< M 



and 



K(s,)K(s2)---K(sn)Gi{x,v,0 



< M 



(59) 



(60) 



for all x,v,n. Then for any arbitrary number of steps in the calculation of 
( 1571) one would obtain a finite value if 



8 max 



A 



(i) 

x,v,t 



rrii 



M < 1 



(61) 



In conclusion: 



Theorem 3 If the smoothness and boundedness conditions f5^) - [61\} are ful- 
filled, the tree-indexed stochastic process Z ^x, f ,tj yields a stochastic solu- 
tion of the configuration space Poisson- Vlasov equation in an external mag- 
netic field. 
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3 Remarks and conclusions 



1) The stochastic solution results established for the Fourier-transformed and 
the configuration space Poisson-Vlasov equations in an external magnetic 
field may, as discussed in the introduction, provide adequate algorithms for 
the parallel computation of localized solutions. That implementation of such 
algorithms is feasible has been shown in Ref.|19] for the stochastic solutions 
associated to branching processes and multiplicative functionals. For the 
Fourier-transformed solutions developed here the algorithms would be quite 
similar, the main difference being the slightly more complex exponential pro- 
cess. However, this extra complexity pays off in allowing for solutions without 
upper time bounds. 

For the tree-indexed processes, that construct the configuration space so- 
lutions, the implementation could lead to larger computer time requirements 
because, for each realization, one has to compute the iteration in Eq. (!57|) and 
then to average over many realizations. 

2) In plasma phenomena in strong magnetic fields there is a hierarchy of 
well separated time scales, the Larmor time scale, the bounce time scale and 
the drift time scale. Separation of the Larmor time scale led to a beautiful 
body of theory that goes by the name of gyrokinetics|28]. A practical mo- 
tivation for the gyrokinetics reduction comes from the possibility to reduce 
the dimension of the numerical codes from 6 to 5 or 4 dimensions. With 
the present improvement of multiprocessor computer power this motivation 
has somehow become weaker, especially because of the additional complexity 
of the gyrokinetics equations if one wants to go beyond the leading order. 
That, to obtain any reasonable accuracy, higher gyrokinetic orders should be 
included in the numerical calculation, is indeed to be expected in view of the 
fact that the exact invariant associated to the gyrokinetics reduction may, at 
best, be obtained by a Borel-summable infinite series [29] 

Nevertheless, if a reduction of the Larmor time scale is desired, the 
stochastic solution approach developed in this paper might also provide an 
appropriate framework for this reduction. Notice in particular that, in the 
configuration space stochastic solutions, the magnetic field evolution acts 
only on the function arguments, that is, on the labels of the stochastic pro- 
cess not on the process itself. Then, averaging techniques or scalar function 
mappings would provide an alternative formulation of gyrokinetics. 

3) In the stochastic solutions for the configuration space equation and for 
the non-uniform magnetic field case, operator labels associated to the parti- 
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cles generated in the tree are carried over and applied to the initial conditions 
when the process arrives to time zero. This entails some additional complex- 
ity in the calculation of the functionals and in the smoothness requirements 
to be imposed on the initial conditions. The need for these operator labels 
arises from the singular nature of the propagation kernels derivatives. A sim- 
ple (one-dimensional) example illustrates this point. Let us assume that a 
probabilistic interpretation is to be given to an integral containing the factor 
dvfi {x, V, t). Then we may replace it by 

~ j 5' iv-v')f,ix,v',t) dv' 

but it is not possible to absorb 6' {v — v') into a probability kernel unless 
some limiting approximation is used 

f I 1 {v-v'f 

i 2sign (t> — v') \v — v'\ lirn W — -e ^ fi {x, v', t) 

with sign(f — v') in the coupling constant and the rest in the probability 
kernel. However, the computation of the approximation entails numerical 
instabilities and to keep the derivative as an operator label seems to be a 
more robust procedure. 

A completely different situation occurs if the derivative of the propagation 
kernel is smooth. This is the case in the Navier-Stokes equation[T6], where 
by an integration by parts the derivative of the heat kernel is controlled by 
a majorizing kernel and absorbed in the probability measure. 
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